Active matter refers to a variety of systems ranging from, e.g., cytoskeleton of living cells on a microscopic scale to schools of fish on a macroscopic scale.^[1] The common feature of these systems is that the constituent units can move actively (self-propelled) by consuming energy. Such a non-thermal motion is the origin of their out-of-equilibrium nature, which leads to very rich intriguing phenomena, e.g., giant fluctuations,^[2,3] spontaneous phase separation,^[4,5] and pattern formation.^[6,7] In addition to the biological active matters, various artificial active systems have been designed and fabricated, which allow the studies of these nonequilibrium phenomena in well-controlled cases.^[2,8–11]

A boundary is ubiquitous in both natural and laboratory active systems. Recent efforts^[12–21] have been devoted to understanding the distributions and dynamics of active units confined by hard boundaries with various shapes. The boundary can block and rectify the motions of these active units and cause them to accumulate nearby. The accumulation is greatly influenced by the shape or local curvature of the boundary.^[15,18] The non-thermal pressure exerted by the active units in turn can drive the immersed objects with asymmetric boundaries to move directionally.^[22–25] In the biological world, the living organisms usually operate with soft boundaries, e.g. molecular motors and microtubules confined by cell membrane, pathogens migrating through the extracellular matrix of the host. Therefore, understanding the interplay between active unit and soft boundary is of fundamental importance.^[26–28]

In this paper, we address this problem by studying a simple two-dimensional (2D) model, in which self-propelled disks are confined inside a closed deformable boundary, which is mimicked by spring-connected beads. This system is relevant to the cases where active organisms are enclosed by membranes, other biofilms or flexible surfaces such as the liquid-water interface.^[29] It can also be realized by confining microscopic self-propelled particles (SPPs) inside a giant vesicle or using elastomer as the boundary for macroscopic SPPs. We find that persistent motion drives the particles to accumulate near the boundary,^[12,13] and dramatically deform its shape through collisions. On the one hand, the boundary deforms when the inhomogeneous aggregation of the particles takes place. On the other hand, the trend of particles accumulating at the place of large curvature,^[18] reinforces such a deformation. In the present paper, we are concerned with the variations of the size and shape of the soft boundary and its coupling with the distribution of SPPs.

We consider a 2D system where self-propelled particles (treated as disks) are enclosed by a soft boundary. Each particle has mass m, and diameter σ. In addition to the thermal Brownian motion at temperature T, each particle is propelled by a force with the amplitude F through its center. The orientation of the force, rotates with time due to fluctuations. The soft boundary is mimicked by a ring of spring-connected passive beads () which are set to be identical to the active particles except being passive, i.e., without propelling force on them. All non-bonded interactions between particles and beads are given by the isotropic and purely repulsive Weeks–Chandler–Andersen (WCA) potential,^[30] with a cutoff at , beyond which , ε is the interaction strength. The bonded harmonic potential between boundary beads is . We set the spring constant and equilibrium bond length to prevent the active particles from crossing the boundary. In the natural circular shape, the perimeter of the boundary , the radius and the enclosed area .

The motions of particles and beads are described by the Langevin equations as follows:

We use the LAMMPS software to perform the simulations.^[31] A square box of with a periodic condition in both x and y directions is adopted. Reduced units are used by setting , , and . is the corresponding unit time. We set ε = 10 and γ = 10, which is large enough so that the motions of particles and beads are effectively overdamped. We treat as an independent parameter because in many cases the rotational noise is athermal.^[12,15] To highlight the influence of the particle activity, we set to be a small value of , inverse of which quantifies the persistence time of the active motion. For each case, it runs in a minimum time of in time steps of .

We consider three initial area fractions of confined SPPs, ϕ = 0.05 (low), 0.2 (medium), and 0.4 (high), where . Figure 1 shows the variations of the typical instantaneous morphologies with the increase of the propelling force on the particles for ϕ = 0.05 and 0.4. Apparently, in both cases, there are two regimes in the variation of the morphology with increasing propelling force. i) For small propelling force, the distribution of particles inside the closed boundary changes from homogeneous distribution to radially inhomogeneous aggregation, i.e., a concentration hole gradually forms in the center due to enhanced aggregation of particles near the boundary; the shape of the boundary roughly keeps circular; and apparently the local fluctuation of the boundary is increasingly suppressed. ii) For large active force, the distribution of particles becomes both radially and angularly inhomogeneous; the boundary is gradually elongated and deviates from a circle; this global shape deformation is the consequence of the inhomogeneous distribution of the particles along the boundary and the trend of their accumulating around the places of high curvature,^[18] which forms positive feedback in the deformation of the boundary. In the following, we quantitatively analyze the variations of the size and shape of the boundary, the inhomogeneous distribution of particles, and their correlations.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Snapshots of the morphologies for the initial particle area fraction ϕ = 0.05 (left) and 0.4 (right) and varying propelling forces. The red lines represent a soft bead-spring boundary and the green spheres denote active particles.

First, we calculate the mean radius of gyration, , of the boundary, based on the positions of its constituent beads, as a measure of its size for a ring polymer (Fig. 2(a)). In two dimensions, the radius of gyration is defined as , where and are the two eigenvalues of the gyration tensor,

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. (color online) Variations of dimensionless mean radius of gyration of the soft boundary with (a) propelling force and (b) the initial particle area fraction, where R0 is the natural radius of the boundary. The inset in panel (a) shows the magnified parts of curves to show the data at small forces. (c) Linear–log plot of the ratio of the mean radius of gyration of the boundary to the persistence length of a free active particle, versus propelling force which quantifies the degree of confinement.

At small F (inset of Fig. 2(a)), we find a crossover from rapid to slow expansion of the boundary with the propelling force around . The initial rapid expansion is attributed to the aggregation of particles to the boundary, which effectively increases the average number of particles that exert a force on the boundary. We calculate the probability distributions of for ϕ = 0.05 under different forces (Fig. 3(a)). The curve of F = 1 shows the narrowest distribution indicating the strong suppression of the local fluctuation of the boundary. This is further proved by the calculation of the relative fluctuation,

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) (a) Probability distributions of the radius of gyration at various forces for ϕ = 0.05, and (b) variations of relative fluctuations of the radius of gyration with propelling force.

To quantify the global shape deformation of the boundary, we calculate the asphericity (deviation from a circle)^[32,33] in Fig. 4. The asphericity in two dimensions is defined as

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) (a) Variations of mean asphericity of the boundary with active force; the inset shows the data at small forces. The decay of the time correlation function of asphericity for concentrations ϕ = 0.05 (b) and 0.4 (c).

Time correlation function of asphericity is a way to quantify the varying rate of the global shape (Figs. 4(b) and 4(c)). We find that the force-dependence of the decay is different for the low and high initial particle concentrations. For ϕ = 0.05, the decays faster (i.e., the shape varies faster) at a larger force (Fig. 4(b)). The change of the decaying rate is appreciable. Taking the time at which as the characteristic correlation time, we have that the correlation time changes from for F = 5 to for F = 80. On the contrary, the decay of becomes slower upon the increase of force for ϕ = 0.4 (Fig. 4(c)), but apparently the change is small. The correlation time changes from for F = 5 to for F = 80. There are two competing factors in determining the rate of the shape variation: expansion of the boundary and the activity of the particles, both of which are enhanced with the increase of active force. The large size of the boundary slows down the shape variation while the large activity of particles leads to the rapid relocation and hence the fast shape variation. The expansion of the boundary with active force is weak for ϕ = 0.05 (Fig. 2(a)) and the enhanced activity of particles dominates and leads to faster decay of at a larger driving force. In contrast, the expansion of the boundary is significant for ϕ = 0.4 and causes a much slower variation of the shape at a larger force. Notice that the correlation time for ϕ = 0.05 is larger than that for ϕ = 0.4 when . We examine the trajectories and find that the possible reason is the different mechanisms for the particle relocation. For ϕ = 0.05 the layer of particles at the boundary is thin and the particles are relocated mostly through the motion along the boundary, while the layer of particles at the boundary is thick for ϕ = 0.4 and besides the motion along the boundary many particles near the center move across the central hole, which can definitely speed up the change of the angular distribution of particles.

To investigate the radial and angular inhomogeneous distributions of the enclosed particles and their coupling with the deformation of the boundary, we calculate the Gini coefficients^[12] and

Fig. 5.

Figure Option ViewDownloadNew Window

Fig. 5. (color online) (a) Variations of radial Gini coefficient, , with propelling force. The inset shows the schematic diagram of partitioning the space radially (left) and angularly along the boundary (right) for calculating the radial and angular Gini coefficients, respectively. (b) Variations of angular Gini coefficient, , with propelling force. The inset shows the magnified parts of curves at small forces. (The standard deviations of the data in panels (a) and (b) are less than 0.09 and the standard errors of the mean are less than .)

Finally, we calculate the local curvature κ of the boundary and also the local pressure on it by statistics. The local curvature , where () is the coordinate of a boundary point, and and are the second and first derivatives with respect to x. The local curvatures are calculated at the bead positions along the boundary in the discrete form of the derivatives. The local pressure exerted by active particles is evaluated by the local force along the local normal direction of the boundary divided by the local arc length. The statistical data for ϕ = 0.05 are shown in Fig. 6(a). It is found that the local pressure is almost independent of the local curvature at small forces (), which corresponds to the small asphericity of the boundary as a whole (Fig. 4(a)). For , local pressure increases with the increase of the local curvature, which accords with the observation on fixed boundary.^[15,18] The correlation between local pressure and local curvature is directly visualized in Fig. 6(b), which shows the heat map of the local pressure and curvature along the contour of the boundary for the case of F = 80 and ϕ = 0.05.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) (a) Variations of average pressure with local curvature κ for ϕ = 0.05. (b) Profile of the instantaneous (i.e., the data are extracted from a typical snapshot) local curvature (inside loop) and local pressure (outside loop) along the boundary for the case of F = 80, ϕ = 0.05.

Boundaries are ubiquitous in both natural and laboratory active systems. Understanding the interplay between active agent and boundary is important for both basic science and applications. Although the living organisms usually operate with soft boundaries, the properties of active agents confined by flexible boundaries have been less addressed. In this paper, we study a very simplified 2D model of SPPs confined by a ring-polymer-like boundary. The cooperative phenomenon between boundary fluctuation/shape deformation and particle distribution/collective motion is found. In summary we find that i) the average size of the boundary varies roughly linearly with the propelling force and the number of enclosed particles; ii) there are two typical regimes in the variation of the morphology with propelling force; iii) there is strong coupling between the shape deformation of the boundary and the particle inhomogeneous distribution; iv) the crossover from regime one to regime two is abrupt for small ϕ but is gradual for large ϕ. The mechanism is similar to the case studied by Nikola et al., who found that the pressure inhomogeneity causes a modulation instability of a short semi-flexible filament.^[34]

Our model is easily realizable in a granular experiment. Also, in a biological experiment, it shows that colloidal particles can be confined inside a fibroblast cell through phagocytosis which suggests a possible way of realizing our model in a biological system.^[35] As a kind of active colloidal cell,^[36] this integrated system has the potential to perform some functions as a machine. Our results provide some insights into the control of the shape by, for example, the amplitude of propelling force and particle concentration, which may help design and manipulate the colloidal machine. Moreover, our model captures some essential features of a biological cell: active motion and deformability. It will be interesting to explore the emergent properties of the assembly of these active colloidal cells.

The very first version of this work was uploaded to the arXiv.org in 2015 and cited by the work,^[37] in which the authors studied the shape and displacement fluctuations of a soft vesicle filled with anisotropic Run-and-Tumble particles.